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We introduce a model description of femtosecond laser induced desorption at surfaces. The sub¬ 
strate part of the system is taken into account as a (possibly semi-infinite) linear chain. Here, being 
especially interested to the early stages of dissociation, we consider a finite-size implementation of 
the model (i.e. a finite substrate), for which an exact numerical solution is possible. By time-evolving 
the many-body wavefunction, and also using results from a novel time-dependent DFT description 
for electron-nuclear systems, we analyse the competition between several surface-response mecha¬ 
nisms and electronic correlations in the transient and longer time dynamics under the influence of 
dipole-coupled fields. Our model allows us to explore how coherent multiple-pulse protocols can 
impact desorption in a variety of prototypical experiments. 


PACS numbers: 78.47.J-, 78.20.Bh, 79.20.La, 31.15.ee 


I. INTRODUCTION 


Femtosecond (fs) laser technology has revolutionised 
our understanding and control of chemical reaction^. 
However, it also presents serious theoretical challenges, 
since ultrafast measurements probe atomic/molecular 
scales far away from equilibrium, where an accurate de¬ 
scription of the concerted motion of electrons and nuclei 
is indispensable to interpret the e^^riment. Progress 
has been made for free molecules^J®, but a compara¬ 
ble understanding is still lacking for surfaceJ^. This is 
unfortunate, since many important reactions are catal¬ 
ysed by a surfac^^Hm, e.g. photocatalytic processes di¬ 
rectly involving light-matter interact iorP^. A ke y tech- 
nology for ultrafast studies are Ti: Sapphire laserffE^E^ 
With a central frequency of SOOnm and pulse durations 
from hundreds to a few fs, these lasers have been piv¬ 
otal to the emergence of no vel su rface-sensitive ultrafast 
microscopies/ spectroscopie^^^^ and to new results on 
chemical reactions and desorption.l^^MUl^ 

An accurate first principle description of ultrafast dy¬ 
namics at surfaces is within the scope of comprehensive 
approaches such as the n on-equilibrium Green’s func¬ 
tion (NEGF) methocP^^ and Time-Dependent Density 
Functional Theory (TDDFT)^®, where there is ongo¬ 
ing effort in this direct iorP^^. However, with current 
treatments of electron-electron and electron-nuclear in¬ 
teractions, the inherently non-perturbative situation of 
desorption is in general not adequately described even 
in the initial stages. Thus, it remains highly relevant 
to theoretically explore simplified models which simulate 
experiments with fs lasers, and utilize the p ossibilities of 
control offered by their pulse structnrd ^^ l ^^ i 

Motivated by this, we present here a model approach 
to pump-probe real-time dynamics of adsorbates, incor¬ 
porating electron interactions, core-hole relaxation, plas- 
mon screening and anharmonic nuclear dynamics. Our 
investigation focusses on the early stages of dissociation 
dynamics; however, in the rest of the paper, this spe¬ 


cific sub-regime we be simply referred to as ’’desorp¬ 
tion” . Our novel description merges elements from three 
popular surface-physics models: the Anderson-Newns- 
Grimley model of chemisorptiorP^M^, the charge-transfer 
model of core photoemissiorP^®^ and the Shin-Metiu 
model of nuclear motiorP^. By considering finite systems, 
we exactly and simultaneously address several competing 
time-scales and response mechanisms, to gain robust, al¬ 
beit qualitative, insight in the ultrafast regime. 


Our main results are: i) For short pulses, it appears a 
unified treatment of electrons and (light) nuclei is needed 
already for the early stage of the dynamics, ii) Desorp¬ 
tion can be controlled in experimentally viable pulse pro¬ 
tocols by manipulating the plasmon response, iii) A mul¬ 
ticomponent TDDFT description of adsorbate dynamics 
unveils highly non-trivial features in the Kohn-Sham po¬ 
tentials. Overall, the results show that our model pro¬ 
vides insight into a wide range of situations for adsorbate 
dynamics (accessible by measuring e.g. the adsorbate- 
surface bond length or the change in desorption yield), 
and represents a novel and versatile benchmark for more 
realistic theoretical treatments. 


The paper is structured as follows: In Section [H] we 
define the Hamiltonian of our model, and discuss its pa¬ 
rameters and the limitations of its finite size realization. 
In Section [TTT| we present results pertaining to the equilib¬ 
rium properties of the system, such as energy level struc¬ 
ture and spectral and response functions. In Section |W 


we discuss general dynamical properties of the system, 
and show how different pulse protocols can be utilized 
to control the evolution of the adsorbate. To get insight 
into the results we in Section [V| present a novel version of 
TDDFT and show the exact Kohn-Sham potentials for 
both the electrons and the nuclei. 
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FIG. 1: The model of Eq. 0’ for a mobile adsorbate (with 
one core and two interacting valence levels) on a 5-site sub¬ 
strate. The adsorbate charge fluctuations are coupled to a lo¬ 
cal/surface plasmon. The substrate-adsorbate effective poten¬ 
tial and the laser perturbation are also schematically shown. 


II. THE MODEL 

The system we consider is illustrated in Fig. and 
consists of a rigid chain (the substrate) of L sites with 
one orbital per site and with a mobile adsorbate at one 
end. The Hamiltonian is given by 

H{t)=H,+Ha + Has + Ht), (1) 

with Hs describing the substrate, Ha the adsorbate, Has 
the adsorbate-substrate interaction, and where A(t) is the 
external laser field. To substrate Hamiltonian is taken to 
be 


Hs = -V ^ ^ujpb'^b ( 2 ) 

{RR'),a 

where c\^ ^ creates an electron with spin a at site R of the 
substrate (we use the index S for the “surface” site of the 
substrate), V is the hopping amplitude in the chain and 
b^ is the creation operator of a plasmon with frequency 
Up. We consider only nearest neighbor hopping, which is 
indicated by the braces around the summation indexes. 

The adsorbate Hamiltonian Ha is given by 

Ha — T 2^^ / j ^v'^v^cT’) 

V,(J 

T ^ ^ UyylTly^a^V' ,(j' W {\ flo^Na (3) 

vv' ,crcr' 

where x and p denote the position and momentum op¬ 
erators for an adsorbate with mass M, and the opera¬ 
tor aj ^ creates an electron with spin a and energy 
in the valence orbital v of the adsorbate. We denote 
by n„,o- = fiR^a = and he single level 

number operators, the latter for a structureless core level 
of energy Cc on the adsorbate, and introduce the total- 
number operator = Ac + a of the adsorbate 

and its ground state value {Na)o' The adsorbate has 
two valence levels (with one exception discussed in re¬ 
lation to Fig. where electrons interact mutually with 
strength (with Uu = 2Ui2 = U 22 = U). In case 
of core-hole photoemission or Auger recombination, the 
valence electrons experience an additional interaction w 


(acting as a local potential), which depends on the core¬ 
level occupatioiP^. Since the system is always in a state 
with 0 or 1 core electronwe have H = H{nc), and 
we consider the case of A^e = A + 1 spin-compensated 
electrons in the other orbitals. 

The adsorbate-substrate interaction Hamiltonian Has 
is given by 

Has = ^ + h.C.) 

v,a 

+ ^(^Na-{Na)o){b^ + b) (4) 

where the first term gives a repulsive ion-ion interaction, 
and the second term is attractiv^^ and due to electron 
hopping between the adsorbate and surface sites, whose 
probability decays exponentially with distance^. To¬ 
gether they create a Morse-like potential landscape for 
the adsorbate, and we consider the parameters r, X and 
g as phenomenological to give a reasonable binding en¬ 
ergy El), vibrational frequency Uph and effective hopping 
amplitude Ve = 

The effective plasmon response to charge fluctuations 
on the adsorbate is given by the last term in Has, where 
the value of {Na)o is to be found self-consistently and 7 
determines the electron-plasmon coupling strength. 

The external field A(t) depends on the experiment con¬ 
sidered, but is restricted by A = 0 for t < 0 . In the dipole 
approximation 

A(t) = ^ ^ Ayy' (t^CLy^Cly'fj , (b) 

Vyi^v' ,a 

and induces transitions between the valence levels. Core¬ 
hole photoemission is treated in the sudden limit (see 
e g| 58 | 59 |j^ where A(t) moves the core occupation from 1 
to 0 at time Tc, to mimic the promotion of a core electron 
to a continuum state (left out of the explicit description). 

Parameters used in the model can be obtained for spe¬ 
cific systems using first principle simulations or experi¬ 
mental observations. As a result, a possible use of our 
model is to address qualitative features in specific realis¬ 
tic systems and experimental conditions. Here, however, 
we are interested in demonstrating the generic usefulness 
of the model for addressing the early stages of surface- 
adsorbate dynamics, and have therefore chosen typical 
physical parameter values, as discussed in the following 
subsection. 


A. Parameters 

To set our energy scale we in the following take V = 1, 
whose value is typically H 1 — 2 eV in a metaP^^. 
The plasmon frequency can vary substantially, but typ¬ 
ical values as measured experimentally are in the range 
4 — 10 eV for surface mode^^, consistent with the clas¬ 
sical results Up = ^irne^/m and Us = Upl\f2 for bulk 
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and surface modes respectively. In the following we take 
ujp = A corresponding to a typical frequency. 

Typical values of the diagonal interaction integrals 
Uyy = U a re in the range 1 — 10 eV for adsorbed 
atom^^^E^EU^ with the off-diagonal elements slightly 
smaller. Increasing the interaction strength enhances the 
degree of correlation among the electrons, and a com¬ 
mon indicator of strong correlations is when the ratio 
U/W becomes larger than unity, with W the bandwidth 
of the system (for the one-dimensional Hubbard model 
W = 4t). In this paper we consider U = 1 and 1/ = 4, 
where the latter gives U/W = 1, i.e it marks the onset 
of the strong correlation regime. 

The energy levels Cy could in principle be taken both 
below, both above or one below and one above the Fermi 
level of the chain, depending on the system of interest 
(corresponding to anionic, cationic or neutral adsorption 
respectively^, assuming they would be filled for an iso¬ 
lated atom). Here we take as a possible (and plausible) 
value ey = — 1//2 for the lower level and adjust of the 
upper level to obtain half-filling on the adsorbate in the 
ground state. The core level is assumed to be deep, and 
its role comes from the Coulomb stabilization energy w 
that acts as a local potential on the valence levels, and 
which is typically 5 —10 eV in magnitud^^. Here we have 
taken w = 6 in order for core level emission to have sub¬ 
stantial impact, as is the case in many naturally occurring 
situations. The role of the mass M is to set the time- 
scale of the nuclear dynamics, and should in principle be 
adjusted depending on the species of atoms considered. 
Since we focus here on qualitative features we have cho¬ 
sen the rather small value M = 352, roughly corresponds 
to the mass of Hydrogen, in order to have reasonable 
simulation times and in line with earlier approacheJ^. 

Our simulations were performed with the bare values 

= 0.3, = 6 and A' = 2 , after which the interatomic 

coordinate is rescaled in order to measure length in units 
of the equilibrium distance Xeq. For a chain with L = 5 
and U = 4, this is equivalent to using the values n = 
^— 2.42, g = = 1.83 and A = X'x^q = 

1.19, the correspond to a binding energy E^j 1.5 (or 
equivalently Eij 1.5 — 3 eV), a phonon frequency ujph — 
0.24 and an effective adsorbate-surface hopping Ve — 1.8, 
that are values typical for chemisorption in the surface 
molecule limit. We take 7 = 1 to have an image potential 
shift of 1 ; = = 0.25, corresponding to intermediate 

coupling regime. 


B. Limitations of the model 

Our model is subject to some limitations, the most 
obvious one being i) the finite size of the substrate, con¬ 
taining only a limited number of de-excitation channels 
(this however, besides making possible an exact solution, 
can have direct relevance for dynamics on thin filmJ^. 
Additionally, ii) Auger recombinatiorP®^ is not consid¬ 
ered. iii) Surface-adsorbate hopping induced by a sur- 
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FIG. 2 : Adiabatic potential energy surfaces (PES) for the 
model system of Eq. <0 with one core and two valence levels 
on the adsorbate. The heat map is for a system with L = 5, 
1/ = 4, ei = — 2 and 62 = 1, with colours related to the den¬ 
sity of PES ranging from high (red) to low (blue). The light 
green curves show the eight lowest PES, with nuclear energy 
levels for the ground state potential surface and a parabolic fit 
(black curve) superimposed. Eor comparison, the lowest PES 
for L = 1 are also shown, for adsorbate interaction U — 1 
(black dashed curves) and F = 4 (black solid curves). Eor 
1/ = 1, the valence levels are at ei = —0.75 and 62 = 2 to 
provide adsorbate level fillings similar to U — A. 


face plasmon should also be included (see e.g.lini) , es¬ 
pecially for ionic chemisorption, iv) The plasmon cou¬ 
pling 7 should depend on the adsorbate distance x, intro¬ 
ducing an electron-nuclear-plasmon coupling term in H. 
v) The electronic interactions are kept local, but longer 
range interactions between substrate and adsorbate can 
in fact play an important role, vi) Lattice vibrations and 
electronic interactions in the substrate are not included. 
Avoiding ii-vi) corresponds to easy but computationally 
demanding extensions and is deferred to future work. On 
the other hand, to avoid i), a semi-infinite substrate can 
be taken into account via e.g. NEGF or TDDFTI^. How¬ 
ever, in this case approximate treatments of interactions 
usually need to be introduced. The finite-size version of 
our model then provides a natural exact benchmark to 
such treatments. 

Finally, typical of real-time dynamics approaches (and 
even for semi-infinite substrates) with finite-timespan 
simulations it is not possible to fully exclude that at long 
times an atom re-adsorbs onto the surface after substan¬ 
tial bond stretching. However, the displacement of the 
adsorbate in the sub-picosecond regime (as investigated 
here) is a prerequisite for complete desorption at later 
times, and the qualitative trends seen in this phase should 
also be reflected in desorption measurement^^^^^. 
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FIG. 3: The imaginary part of the local Green’s function Gw 
and density-density response function Xw is shown for the 
two valence levels vi (lower) and V 2 (upper). All panels show 
the same quantity without (shaded regions) and with (lines) a 
local plasmon mode of frequency cjp = 4 and coupling 7 = 1 . 
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FIG. 4: The imaginary part of the local Green’s function 
Gw is shown for the two valence levels vi (lower) and V 2 
(upper). The shaded regions show the case of an inhnitely 
heavy adsorbate and the lines correspond to a mass M = 352, 
for a plasmon mode of frequency cjp = 4 and coupling 7 = 2 . 


nine and are well separated, for the interaction strengths 
U = A (solid lines) and U = 1 (dashed lines). The dis¬ 
sociation energy (the difference between minimum and 
asymptotic values of the lowest PES) is lower in the for¬ 
mer case, which is also true for L = 5. For the parameter 
regime we consider, the plasmon has minor impact on the 
PES structure, while crucially affecting the short-time 
dynamics. To gain further insight into the role played 
by the plasmon, we in the following section discuss the 
spectral and response functions. 


III. THE EQUILIBRIUM CASE 

In this section we discuss some equilibrium properties 
of the system, for the parameters L = 1 and 5 and U = 1 
and 4. To find the exact ground (and initial) state \g) 
we diagonalise il(0) in the basis ^/c, ^ 5 )}, where the 

n^cris are site/orbital/spin occupations, Xk denotes the k- 
th mesh point on a uniform grid, and is the plasmon 
occupation number. 


A. Potential energy surfaces 

The heat map in Eig. [^represents the density of poten¬ 
tial energy surfaces (PES) for L = 5, obtained via bin¬ 
ning around each adsorbate position (for each Xk^ there 
are 3775 PES). We see a large number of surface crossings 
starting already at low energies, that begin to merge at 
larger energy to form a quasi-continuum. The black curve 
shows the harmonic approximation to the ground state 
PES (corresponding to a phonon energy ujph = 0.24) 
that however breaks down almost immediately, as it can 
be seen from the difference in the lowest nuclear energy 
levels of the real and harmonic PES (black and green 
horizontal lines respectively). Eor comparison we show 
results for a dimer, where the number of PES are only 


B. Spectral and response functions 

We here expand on the importance of plasmon effects 
in equilibrium, and to extract this information we look 
at the spectral functions A{uj) and B{uj) of the local one- 
particle Green function G and density-density response 
function x for the adsorbate levels. The zero-temperature 
Green function in the site basis is defined in equilibrium 
as 

Gij{t) = I(V’|r{ci(t)c](0)}|V’) (6) 

where T is the time-ordering operator. Taking the 
Eourier transform of this expression we find the spec¬ 
tral function to be A{uj) = 2isgn{uj)^G{uj). Similarly 
the density-density response function is defined by 

Xij{t) = j{ip\T{Ahi{t)Ahj{0)}\ip), (7) 

where Ahi = hi — rii is the local density fluctuation op¬ 
erator, and from which the spectral function B{uj) = 
2isgn(co’)^x(co’) can once again be found by Eourier trans¬ 
form. 

The electron-plasmon interaction term can be re¬ 
moved from the Hamiltonian via a Lang-Eirsov 
transformatioii^, with the effect of renormalizing the 
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values of the onsite energy and the electron-electron in¬ 
teraction. More specifically, we expect shifts of the sort 

Uyy> “ 2^^ jU^ . ThCSC 
features can both be observed in Fig. above, where 
the interaction induced gap in is diminished, and the 
single-particle levels in are shifted towards lower en¬ 
ergies. 

For 7 = 2, ^G shows several distinct peaked features 
that, loosely speaking, can be s een as emerging from 
broadened plasmon satellite^^^^, where the broadening 
is largely due to the mobility of the adsorbate and the in¬ 
herent fluctuations. To support this argument, in Fig. 
we compare the density of states for an infinitely heavy 
adsorbate to a system with M = 352. As soon as the 
mass becomes finite, i.e. the adsorbate mobility and the 
position fluctuations are increased, each peak is broad¬ 
ened and split into several smaller ones, and the weight 
of the distribution is shifted towards higher energies. 


IV. DESORPTION DYNAMICS 

Starting from \g), for t > 0 the exact many-body wave- 
function is time evolved via the short iterated Lanczos 
algor it hnP^. To induce desorption dynamics we apply 
single- or double-pulse fields, and look at the electron 
density n^(t) at the adsorbate level u, the mean internu- 
clear position x{t) and the nuclear probability distribu¬ 
tion Pt{xk)- The pulses we consider have a FWHM of 6fs 
(30fs), and a carrier wavelength of 800nm corresponding 
to a photon energy uj eV. 



A. Dependence on size and interaction strength 

To assess the role of the substrate on the adsorbate 
dynamics, in Fig. we show the adsorbate wavepacket 
Pt{xk) for substrates of length L = 1, 3 and 5 and van¬ 
ishing electron-plasmon coupling (red, orange and green 
curves). After the pulse has been applied, Pt{xk) pro¬ 
gressively spreads over larger internuclear distances. For 
L = 1 the wavepacket is clearly split into two sepa¬ 
rate structures, while for larger L it is more uniformly 
distributed. We define a desorption yield according to 
Y{t) = 1 — Pt(x)dx, where xq is chosen as the small¬ 
est value for which T(0) = 0. At large times Y = 0.65, 
0.85 and 0.63 for L = 1,3 and 5 respectively, and a 
non-monotonic Y appears to be a rather general fea¬ 
ture; we also observed it for much longer, non-interacting 
chains in the Ehrenfest approximation. In the dipole ap¬ 
proximation Y appears to be only mildly sensitive to L. 
This is partly due to the system approaching the sur¬ 
face molecule limit 1.8 > V), and partly 

because the excitations induced by A are within the ad¬ 
sorbate: we observe a stronger dependence of Y on L for 
fields coupling to the valence level density (not shown 
here). This suggest desorption scenarios where non-local 
effects from the substrate play only a small role. For 


FIG. 5: Dependence of the adsorbate wavepacket dynamics 
on the length L of the substrate and the interaction strength 
U. Panel a) shows snapshots of the wavepacket for L — 1 
(red) and L = 3 (orange), as well as for L = 5 with an 
electron-plasmon coupling 7 = 0 (green) and 7 = 1 (blue). 
In h) we show the evolution of the nuclear wavepacket for 
= 1 (green) and 1/ = 4 (blue). In both panels A^^/(t) = 
/'^cos(cjt), with to = 10, r = 13 and uj — Ott/S, of 
amplitude A = 2. 


L = 5, we also include the local plasmon. The behaviour 
with and without plasmons (blue and green curves re¬ 
spectively) is about the same; however, with the plas¬ 
mon, Y = 0.59 indicative of less desorption within the 
considered time-interval. 

Electronic interactions are expected to have a non- 
negligible effect on desorption. To illustrate their role 
(Fig. 1 ) 3 ), we compare the distributions Pt{xk) for U = 1 
and 4 that give the respective desorption yields Y = 0.42 
and 0.59. Eor U = 1 the adsorbate electronic density 
(not shown) fluctuates more in time, while for large in¬ 
teraction these oscillations are quenched. Thus the sup¬ 
pression of charge fluctuations to/from the substrate by 
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FIG. 6: Heat map of the occupation in the hundred first nat¬ 
ural orbitals of the reduced nuclear density matrix, as a func¬ 
tion of time. Superimposed is the entanglement entropy Sn 
for the nuclear subsystem. 


electronic correlations appears to affect the probability 
of desorption. 


B. Natural nuclear orbitals 

Using an adequate space grid can be a computational 
bottleneck for calculations addressing desorption; fur¬ 
thermore, differently from on-resonance experiments, ul- 
trashort pulses involve a large spectrum of frequencies, 
and many PES are simultaneously involved. But how 
many is ’’many”? As an empirical answer, the heatmap 
in Figj^ shows the occupation of the natural orbitals of 
the reduced nuclear density operator T{t) = Tre,p/p(t), 
after the application of a 6fs pulse of amplitude A = 3. 
The other system parameters are the same as for the blue 
curves in Fig. [^, Here the trace is taken over the elec¬ 
tronic and plasmonic degrees of freedom. According to 
these results, a calculations would only require the first 
50 or so natural orbitals to keep significant precision, re¬ 
ducing by a factor of 20 the currently used space grid 
basis, but still maintaining full correlation between nu¬ 
clear and electronic degrees of freedom. The entity of 
these correlations can be easily realized when looking at 
the nuclear entanglement entropy Sn = TrTlogT, which 
grows quickly already in the very early stages of desorp¬ 
tion. 


C. Manipulating the system via pulse control 

We now go on to discuss the dynamics of a surface- 
adsorbate system induced by different pulse protocols, 
and to highlight the different time-scales at play we start 
by discussing a very simple experiment where either the 
pulse duration or amplitude is varied. For this pur¬ 
pose we find it convenient to analyze the behavior of the 
the bond kinetic energy Kas = + c|c^), which 




FIG. 7: Time evolution (in units of fs) of the nuclear den¬ 
sity (red/blue), average position (orange/green) and plasmon 
density (back), for different pulse protocols (front). The pa¬ 
rameters are as discussed in Section lira with L = 5 and 
U = 4. In a) two pulses of identical integrated intensity and 
duration 6 fs (orange) and 35 fs (green) are shown; in b) two 
pulses of identical amplitude A = 2 and duration 6 fs (orange) 
and 35 fs (green) are shown. The panels at the bottom show 
the respective bond kinetic energies Kad- 


we take as our measure of the adsorbate-surface bond 
strength. 

In Fig. [Tfi we compare two 800nm (IR) pulses with 
FWHMs of 6fs and 35fs, and of equal integrated intensity. 
In the first case we see a significant amount of plasmon 
excitation and reduction of the bond kinetic energy be¬ 
tween the surface and the adsorbate. In the second case 
there is much less change in both plasmon density and ki- 
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FIG. 8: Time evolution of the nuclear density (red/blue), 
average position (orange/green) and plasmon density (back), 
for two different pulse protocols (front). Parameters are as 
in Fig. In a) two pulses are applied consecutively with 
delays of 14fs or 18fs, and A^^/(t) = cos{ujt) + 

cos{ujt-\-ip). Here to = 10, r = 13, uj — 67r/8 and 
A = 2, with ti = 24 (red/green) or ti = 28 (blue/orange), and 
(f chosen to give the second pulse the same envelope-carrier 
relation as the first. In b) a core electron is removed during 
the action of a 6fs pulse, at Tc = 10 (red/green) or Tc = 11.3 
(blue/orange) with re = 6. 


netic energy, leading also to much less desorption. Since 
the only difference is the duration of the pulse the differ¬ 
ent outcomes are most likely due to the fact that, for the 
shorter pulse, non-adiabatic effects in the response play 
a greater role. 

In Fig.[^ we instead compare two pulses of equal am¬ 


plitude A = 2, but with with FWHMs of 6fs and 35fs. 
The latter leads to a much larger displacement of the ad¬ 
sorbate within the same time-frame, as expected from its 
greater (by a factor 2.5) deposited energy. For the 35fs 
pulse the evolution of the system is dominated by the 
shape of the field, while for the 6fs pulse its behaviour 
is to a larger extent determined by the internal dynam¬ 
ics. While this is a very reasonable physical result, it 
also points to the observation that ultrafast dynamics 
within the first 50 fs occurs both for the electron and 
nuclear dynamics in a correlated fashion, and thus sep¬ 
arating these two time scales might not be appropriate. 
This kind of correlation is experimentally observed for 
molecular systemJ^. 

Finally, we use the model to gain insight into two 
prototypical experiments, based on pulse durations 
and wavelengths realizable with existing state-of-the- 
art laser^^^HzSl Fig. we apply in a pump-probe 
manne]®^ two 6fs pulses with delays differing by 1.5 IR 
cycles, to explore the possibility of coherent control of 
the adsorbate motion. The effect of changing the delay 
is clear: while both cases give an increase in the average 
adsorbate position, the shorter delay leads to a 15 per¬ 
cent larger stretching (after lOOfs). This is due to the 
plasmon: in contrast to the equilibrium case, it acts as 
a strong harmonic perturbation which, during its cycle, 
can be rei nforced by applying a second pulse at the right 
pointUUIll. 

As a second example, in Fig.[^ an instantaneous core¬ 
level photoemission (PE) is combined with an IR pulse, 
as conducted in IR+XUV experiments using high har¬ 
monic generation technology^. Ejecting an electron has 
a huge impact, and the average position of the adsorbate 
in time is almost four times larger compared to only an 
IR pulse (see Eig. IT)- Interestingly, a change of the PE 
time with less than an IR cycle will significantly influence 
both the adsorbate dynamics and the plasmon behaviour. 
As an overall, final observation about Eig. it can also 
be seen that while the moderate field strengths used here 
do not lead to desorption within the fi rst tenths of fs (as 
found in some molecular system^^^^^, significant bond¬ 
stretching occurs so rapidly that separating nuclear and 
electronic timescales may not be justified for these exper¬ 
iments. 


V. A TDDFT PERSPECTIVE 

As a way to obtain insight into the desorption 
dynamic s we consider a multi-component TDDET 
approackP^l^, specialised to electrons on a lattice. A 
system with electrons on a set of spin-orbitals in 

a lattice, and nuclei at {Rm} = (where 

and Cm denote space and spin variables respectively), can 
be described by a Hamiltonian with external potentials 
and e®^^(R, t). On the nuclear side, we choose 
as fundamental variable the diagonal F(R, t) of the one- 
particle density matri:jP^. Eor the electrons we observe 




























FIG. 9: Panel a) gives the exact nuclear potential cks for 
a dimer with one valence level at ei = —Ul2, while b) shows 
the electron density at the adsorbate (red) and the exact elec¬ 
tronic KS potential \Tks\ (blue) and a,Tg{TKs) (green) for the 
same system. The insets are snapshots att = 0t = 17 and 
t = 34) of the nuclear wavepacket and corresponding cks- In 
all cases A^^/(t) = cos{ujt), with to = 10, r = 13 

and Ct; = Gtt/S, of amplitude A = 10. 


that increasing the internuclear distance should result in 
a reduced hopping probability, so the KS Hamiltonian 
Hks should in the lattice basis have (in general complex) 
matrix elements where both modulus and phase 

can vary. This is taken into account by a generalization of 
a lattice time-dependent current DFT that uses the com¬ 
plex bond current as basic variabl^^, defined as Qfj (t) = 
(t )+ Pij{t). For a purely electronic system only 
p^j{t) = {^p{t)\al^ajcr\'ip{t)) would enter, and the explicit 
electron-nuclear coupling K^j({R^}) is contained in the 
second term pF(t) = (' 0 (t)|Kij({Rm}) 4 ^a^v|' 0 (t)) (there 
is of course an implicit depe ndence through the state vec¬ 
tor). Proceeding as a bijective mapping between 

the fundamental variables {Qjj , F) and the potentials 
e®^^) can be establisheci^ and as discussed irP^T- 
representability (and non-interacting T-representability) 
by is ensured when \pij{t)\ > 0 . 

We show in Fig. [^results for a dimer (L = 1 ) without 
plasmons (7 = 0 ), with one valence level on the adsor¬ 


bate, and write R ^ x for the single, ID nuclear co¬ 
ordinate of our model. For the system considered the 
Kohn-Sham Hamiltonians for the electronic and nuclear 
system are respectively 

= E (^f ^ ’ r] (t)c7c,- . + h.C.) (8) 

and 

=T.^ + ^Ks[Qij,T]{x,t). (9) 

k 

and this simplified case (a purely electronic Hubbard 
dimer is non-interacting 'T’-representabl^^ already shows 
essential features of the KS potentials and 

which in fact can be constructed exactly. 

For the electrons, the complex are deter¬ 

mined via numerical reverse engineering. To determine 
^Ks{x,t) we perform an exact factorisation of the wave- 
function of the interacting electron-nuclear systenPH, by 
defining a nuclear wavefunction x(x,t) = t). 

Choosing the gauge where the vector potential is zercP^, 
the exact nuclear potential is e = ( 2 M)“^ [(9^, In^)^ + 
9a,a,ln^ — {dxS)‘^] — dfS. This is also the exact KS nu¬ 
clear potential after introducing a discrete 

mesh Xk for the nuclear coordinate). 

In Fig. it the splitting of the nuclear wavepacket 
during desorption is seen to reflect in the behaviour of 
^Ks{xk^t)\ for t 10 and 17, exs develops dips to dy¬ 
namically push outwards part of the wavepacket. Be¬ 
tween these times it undergoes several dynamical cor¬ 
rugations, related to the electronic oscillations in turn 
induced by the external pulse. This is seen also in the 
snapshots at the bottom of Fig. it, showing the nuclear 
wavepacket and the corresponding exs- The modulus 
T^s electronic KS potential decreases for large 

adsorbate-substrate distances (as expected on physical 
grounds), while the phase appears to grow in a steady 
fashion; as a reference we show the adsorbate electronic 
density (red curve). During the central phase of the 
desorption process we see rapid oscillations in on 

speculative grounds, such a non-trivial temporal pattern 
suggests that, in general, desorption and charge trans¬ 
fer at surfaces can be quite challenging to a TDDFT 
descriptiorP^, e.g. when adiabatic TDDFT treatments 
are considered. 


A. Adiabatic approximation 

The potentials and exs can in general be very 

complicated, but there exist cases where an adiabatic 
approximation could be expected to perform rather well, 
even in the limit of significant desorption. In Fig. p!Q| we 
compare the expectation values of the electronic density, 
bond kinetic energy and mean internuclear position, with 
the external field given by a Gaussian pulse with FWHM 
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FIG. 10: In the three leftmost panels we show the electron density, bond kinetic energy and mean internuclear distance, after 
the application of a square potential with amplitude A = 2 (blue), A = 6 (green) or ^ = 10 (orange). The three central panels 
show the argument and modulus of the exact electronic Kohn-Sham potential as well as the time structure of the 

external laser field. In the panels furthest to the right we show 3D plots of the exact nuclear potential for pulses of amplitude 
A = 2 (top) and A = 10 (bottom). 



Time Time 

FIG. 11: In the three leftmost panels we show the electron density, bond kinetic energy and mean internuclear distance, after 
the application of a square potential with amplitude A = 2 (blue), A = 6 (green) or A = 10 (orange). The three central panels 
show the argument and modulus of the exact electronic Kohn-Sham potential as well as the time structure of the 

external laser field. In the panels furthest to the right we show 3D plots of the exact nuclear potential for pulses of amplitude 
A = 2 (top) and A = 10 (bottom). 


6fs and carrier wavelength of SOOnm, and an amplitude 
of 2 (blue), 6 (green) and 10 (orange) respectively. In all 
cases there are rapid oscillations in the electronic quan¬ 
tities, which are reflected in the KS potentials. Even in 
the case of weak perturbation where the adsorbate stays 
bound and performs oscillations (blue), we see features 
in the potentials that could be hard to reproduce in the 
adiabatic approximation. 


In Fig. [^we show the same quantities for a square¬ 
like external field, given by A(t) = X[o,r](0 sin^ (7rt/2r) + 

X[T,2T]{t) + X[2T,3T]{t) [l - Sin^ {TTt/2T - tt)], where xi{t) 
is the characteristic function of the interval I {xiif) = 1 
for t G / and 0 otherwise) and r = 15. In this case the be¬ 
havior of the KS potentials is much smoother, especially 
during the action of the pulse. For the strongest pertur¬ 
bation, leading to a significant amount of desorption, the 
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Time Time 

FIG. 12: In the three leftmost panels we show the electron density, bond kinetic energy and mean internuclear distance, after 
the application of a square potential with amplitude A = 2 (blue), A = 6 (green) or ^ = 10 (orange). The three central panels 
show the argument and modulus of the exact electronic Kohn-Sham potential as well as the time structure of the 

external laser field. In the panels furthest to the right we show 3D plots of the exact nuclear potential for pulses of amplitude 
A = 2 (top) and A = 10 (bottom). 


amplitude of the electronic potential shows rapid oscil¬ 
lations at the end of the pulse, in contrast to the phase 
which is well behaved for all times. 

In Fig.[^we show the case of a smeared step potential, 
given explicitly by A(t) = X[o,r]{t) {7Ttl2T) + 0{t — r) 
with r = 20. We see that all quantities change smoothly 
during the application of the pulse, independently of the 
strength of the perturbation. For the strongest external 
field we see clear indications of desorption for large times, 
and as expected the amplitude of the hopping parameter 
approaches zero in this limit. It is thus likely that an 
adiabatic approximation within TDDFT would perform 
well under these circumstances. 

As a final comment we note that the nuclear KS po¬ 
tential is typically smoother for stronger perturbations, 
which can be understood from the evolution of the nu¬ 
clear density distribution (see Figs. |5][^. For weak per¬ 
turbations part of the density is ejected around the time 
that the pulse acts, while another part stays bound. The 
part that remains in the potential well performs oscilla¬ 
tions, and each time it reaches the turning point of the 
potential part of the density is emitted. These succes¬ 
sive emissions are reflected in the successive dips in the 
nuclear KS potential, that can be see for all three per¬ 
turbations. For a strong perturbation there is only one 
emission, where the whole wavepacket is released, and 
after this the potential behaves quite smoothly. 

VI. CONCLUSIONS 

In this work, we have introduced an exactly solvable 
model for electron-nuclear dynamics of adsorbates in¬ 


duced by ultrashort laser pulses. Though finite in size, 
the systems we can treat still contain the rich behaviour 
expected for an adsorbate-surface system. To illustrate 
the broad scope of the model, we briefly touched upon 
several issues, e.g. adsorbate dynamics in the surface 
molecule limit, electronic correlations and manipulation 
of plasmon dynamics. We also showed that the model 
can be a valuable aid in devising laser schemes to control 
the outcome of surface studies of light-matter interaction. 
Further, it can be used to benchmark more realistic theo¬ 
retical approaches where approximations are necessarily 
introduced (it was e.g. used here to gain insight into 
general features of the KS potentials of TDDFT during 
desorption). As adsorbate systems are a key paradigm to 
explore light-matter interactions, and new laser sources 
will increasingly be applied to surfaces and thin films, 
work to extend and apply our model in several directions 
is under way. 
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